Giant fluctuations at a granular phase separation threshold 
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We investigate a phase separation instability that occurs in a system of nearly elastically colliding 
hard spheres driven by a thermal wall. If the aspect ratio of the confining box exceeds a thresh- 
old value, granular hydrostatics predict phase separation: the formation of a high-density region 
coexisting with a low-density region along the wall that is opposite to the thermal wall. Event- 
driven molecular dynamic simulations confirm this prediction. The theoretical bifurcation curve 
agrees with the simulations quantitatively well below and well above the threshold. However, in a 
wide region of aspect ratios around the threshold, the system is dominated by fluctuations, and the 
hydrostatic theory breaks down. Two possible scenarios of the origin of the giant fluctuations are 
discussed. 
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I. INTRODUCTION 

Dynamics of a system of inelastically colliding hard 
spheres have attracted a great deal of recent interest 
1, 2], in particular in the context of validity of kinetic 
theory and hydrodynamics of rapid granular flow devel- 
oped in the 80-ies Q • Hydrodynamics looks ideally suit- 
able for a description of large-scale patterns observed in 
rapidgranular flows: a plethora of clustering phenom- 
ena U vortices oscillons |(|, shocks 0, etc., that 
are difficult to understand in the language of individual 
particles. However, a first-principle derivation of a uni- 
versally applicable continuum theory of granular gas is 
not a simple task, even in the dilute limit. The use of 
the Enskog equation, the starting point of a systematic 
derivation of the constitutive relations of granular hydro- 
dynamics, is based on the Molecular Chaos hypothesis. 
This hypothesis is justified for not too large densities and 
for an ensemble of elastic hard spheres. Its use for inelas- 
tic hard spheres is not obvious, as inelasticity of the par- 
ticle collisions introduces inter-particle correlations [g. 
The correlations become stronger as the inelasticity of 
the collisions increases. On the contrary, for nearly elas- 
tic collisions, 1 — r 2 <C 1 (where r is the coefficient of 
normal restitution) the correlations are small, and the 
Enskog equation can be safely used. 

An important additional assumption, made in the pro- 
cess of the derivation of hydrodynamics from the En- 
skog equation, is scale separation. Hydrodynamics de- 
mands that the mean free path of the particles be much 
less than any characteristic length scale, and the mean 
time between two consecutive collisions be much less than 
any characteristic time scale described hydrodynamically. 
This condition should be verified, in every specific sys- 
tem, after the hydrodynamic problem is solved and the 
characteristic length and time scales determined. Again, 
it is safe to say that this condition can be satisfied if the 
particle collisions are nearly elastic |9l ll0llll| . Restrictive 
as it is, the nearly elastic limit is conceptually important 
just because granular hydrodynamics is expected to work 



here. 

Another potentially important, albeit largely unex- 
plored, limitation of the validity of granular hydrody- 
namics (or, rather, of any continuum approach to rapid 
granular flow) is due to the noise caused by the discrete 
nature of particles. Noise is stronger here than in classical 
(molecular) fluids simply because the number of particles 
is much smaller. In addition, noise can be amplified at 
thresholds of hydrodynamic instabilities as found, for ex- 
ample, in Rayleigh-Benard convection of classical fluids 

m 

The validity of hydrodynamic description in general, 
and the accuracy of constitutive relations in particular, 
can be conveniently checked on symmetry-breaking insta- 
bilities that are abundant in rapid granular flows. The 
example of a symmetry-breaking instability that we con- 
sider in this work deals with a very simple setting: a 
two-dimensional (2D) system of nearly elastically col- 
liding hard spheres, confined by a rectangular box and 
driven by a thermal sidewall at zero gravity. The set- 
ting is described in detail in Sec. II. The basic steady 
state here is the "stripe state": a stripe of enhanced 
density at the wall opposite to the driving wall |lfjj |. In 
the continuum language, the stripe state is uniform in 
the lateral direction, by which we mean the direction 
parallel to the driving wall. Within a certain range of 
parameters (delineated below) , steady-state equations of 
granular hydrodynamics predict spontaneous symmetry 
breaking instability of the stripe state, when the aspect 
ratio of the confining box exceeds a certain threshold 
(l3l IT3 . IT5L IT6| . The instability leads to phase separation: 
the development of "droplets" (high-density domains) co- 
existing with "bubbles" (low-density domains). For very 
large aspect ratios of the box, this symmetry-breaking 
instability has been recently observed in event-driven 
molecular dynamic (EMD) simulations, and described by 
a phenomenological continuum model |17|. The present 
work is devoted to a more detailed investigation of the 
phase separation instability in the range of aspect ra- 
tios comparable to the threshold value. We employ, in 
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Sec. Ill, the equations of granular hydrodynamics (or 
rather hydrostatics) to compute the supercritical bifur- 
cation curve for the phase separation instability. Then 
we report, in Section IV, on extensive EMD simulations 
that show that this bifurcation curve is quantitatively ac- 
curate well below and well above the threshold value of 
the aspect ratio. Unexpectedly, the hydrostatic theory 
fails in a relatively wide region of aspect ratios around 
the threshold value, where the system is found to exhibit 
giant fluctuations. In an attempt to get insight into the 
mechanism of this anomaly, we investigate, also in Sec- 
tion IV, the dependence of the magnitude of fluctuations 
on the total number of particles in the system. A sum- 
mary and discussion of our results is presented in Section 
V. 



II. MODEL SYSTEM AND HYDROSTATIC 
EQUATIONS 

Let N hard spheres of diameter d and mass m = 1 
move in a 2D rectangular box L x x L y . The inelasticity 
of particle collisions is parameterized by a constant co- 
efficient of normal restitution r. Particle collisions with 
three of the walls are elastic. The fourth, thermal wall 
is located at i = L x . Upon collision with it, the nor- 
mal component of the particle velocity is drawn from 
a Maxwell distribution with temperature To .lCj, while 
the tangential component of the particle velocity is pre- 
served. 

Working in the nearly elastic limit 1 — r 2 <C 1 and 
employing the Navier- Stokes hydrodynamics we in- 
troduce the number density n(r, t), granular temperature 
T(r,t) and mean-flow velocity v(r,i). Energy input at 
the thermal wall can be balanced by the dissipation due 
to inter-particle collisions. Therefore, we assume that the 
system reaches a zero-mean-flow steady state v = 0, and 
is therefore describable by the simple momentum and en- 
ergy balance equations: 



p = const , V • (k VT) = I . 



(1) 



Here p is the pressure, k is the thermal conductivity and 
/ is the rate of energy loss by collisions. The hydro- 
static Eqs. should be supplemented by constitutive 
relations: p, k and / in terms of n and T. These rela- 
tions are derivable systematically only in the dilute limit 
0,^3- Being interested in moderate densities, we shall 
employ the well-known constitutive relations by Jenk- 
ins and Richman |20|, that account for excluded particle 
volume. In the nearly-elastic limit one can neglect the 
inelasticity correction terms in p and well as the 

small density gradient term, proportional to 1 — r, in the 
heat flux pj . 

Equations Q can be rewritten in terms of a single 
variable: the scaled inverse density z(x, y) = n c /n(x,y), 
where n c — 2/(\/3d 2 ) is the hexagonal close-packing den- 
sity. In scaled coordinates, r/L x — > r, the box dimensions 



become 1 x A. where A = L y /L x is the box aspect ratio. 
We obtain 



V ■{F{z)Vz) = r}Q{z), 
where F(z) = A(z) B{z), 
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and rj = (27r/3)(l — r)(L x /d) 2 is the hydrodynamic in- 
elasticity parameter. Introducing ip (x, y) — J Q F(z')dz', 
we arrive at the following equation: 



(4) 



where Q(ip) = Q [z(tp)] (in the following the symbol ~ is 
omitted). The boundary conditions are 
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Finally, the number of particles is conserved: 
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(6) 



The hydrostatic problem is fully determined by 

three scaled parameters: the area fraction /, 77, and A. 
Notice that the steady-state density distributions are in- 
dependent of To, as the hard sphere model does not in- 
troduce any intrinsic energy scale. 



III. STRIPE STATE, SYMMETRY-BREAKING 
INSTABILITY AND BIFURCATION CURVE 



The trivial steady state of the system is a laterally uni- 
form cluster of particles located at the wall x = 0, oppo- 
site to the thermal wall 153, see Fig. 1. This state will be 
called the stripe state. In the language of hydrodynam- 
ics, it is described by the y-independent solution of Eqs. 
(@J-(|HJ); we shall denote it by z — Z(x), correspondingly 
ip = ^(x). 

It was predicted that, in a wide region of the pa- 
rameter space (/, 77, A), the stripe state should give 
way, by a symmetry-breaking bifurcation (either super- 
critical or subcritical), to a laterally asymmetric state 
[T3L IT1. ITH H(| . For very large aspect ratios A, this 
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FIG. 1: The stripe state for rj = 11,050 and / = 0.025. We 
show the scaled density versus scaled coordinate x obtained 
(a) by solving numerically Eqs. JIJ -© in one dimension 
(line) and (b) in EMD simulation with N = 2 ■ 10 4 parti- 
cles for A = 0.1 (squares). The inset shows a snapshot of 
the system from the EMD simulation (the hot wall is on the 
right). Because of a finite image resolution the particle num- 
ber density in this and other snapshots may look higher than 
it is. 



phase-separation instability has been observed in EMD 
simulations [l^. For a laterally asymmetric steady state 
one can write 



tp(x, y) = *(x) + f n {x) exp(mfcy) , 



(7) 



where ip- n (x) = ip^(x). What happens close to the su- 
percritical bifurcation point? Here the leading terms are 
those with n = ±1, while (po ~ ipf, if2 ~ <p\, <pa ~ 
etc. The bifurcation point itself can be found from the 
linear eigenvalue problem 

Vik -VQ^ Vik - k 2 c ipi k = , (8) 

<p' lk (0) = and Vlk (l) = (9) 
that was analyzed in Refs. [l3L IliL il5j . Here 
Qa,(x) = F^ 1 dQ/dz\ z=z[x) . 

For given rj and /, one obtains the eigenvalue k = k c (ri, f) 
and corresponding eigenfunction tpik(x). The modes 
with k < k c (r), /) are unstable. Within a spinodal in- 
terval fi(i]) < f < f 2(11)1 the effective lateral compress- 
ibility of the gas is negative, and this is the mechanism 
of the instability 0, 0|- At 77 3> 1, there is a range 
of / such that k c and (pik(x) become insensitive to the 
precise form of the boundary conditions at the driving 
wall. This is the universal "localization regime", when 
the eigenfunction ipik (x) is exponentially localized at the 
wall opposite to the driving wall [HI13. The spinodal 
interval exists for r\ c < rj < 00; it shrinks to zero at 



Tj = tj c — 344.3 0, 0| . It has been recently shown, 
for a different boundary condition at the driving wall, 
that the bifurcation from the stripe state to a phase- 
separated state is supercritical within some density in- 
terval /-(??) < / < /+(??), which is located within the 
spinodal interval. On each of the intervals /1 < / < /_ 
and f+ < f < fs, the bifurcation is subcritical |16| . 

As we have already noted, the present work focuses 
on the phase separation via a supercritical bifurcation. 
To obtain the asymptotics of the supercritical bifurca- 
tion curve close to onset, one should go to the second 
order of the perturbation theory and take into account, 
in Eq. (0), the terms n = 0, ±1 and ±2. In this way one 
obtains three linear ordinary differential equations, pre- 
sented in Ref. |16| . where the same problem was solved 
for a different boundary condition at the driving wall. 
The solvability condition for these equations yields 
the bifurcation curve: A versus k^ — k 2 . The amplitude 
A can be uniquely defined by the relation 

<p(x) = A$ (x) + A\A\ 2 Sip(x) , 

where $o(x) is the solution of Eqs. © and 10 such that 
$o(0) = 1, while 5ip(x) = 0(1). This yields 

A(k 2 c -k 2 ) =CA\A\ 2 , 

where C =const. The trivial solution A — describes the 
stripe state, while the nontrivial one, k 2 — k 2 = C\A\ 2 
describes the bifurcated state. The constant C can be 
computed numerically. C > (< 0) corresponds to su- 
percritical (subcritical) bifurcation. We present here the 
resulting bifurcation curve for Y c , the (normalized) y- 
coordinate of the center of mass of the granulate 



Y r = 



Ia dx I-A 2 /2 yn(x,y)dy 
A Jo dx J^l 2 /2 n(x,y)dy 



(10) 



Let us fix r\ and / and treat A as the control parameter. 
When A is slightly larger than A c = 7r/fc c (/), only the 
fundamental mode k = n/A is unstable, and the bifur- 
cation curve has the form 



\Y C \ = T (A - A c ) 



1/2 



(11) 



Here 



T = 



2 3 / 2 /o 



Jo = 2 



dx 



01 



Z 2 F 



and $01 (x) is the solution of initial- value problem for Eq. 
© with the initial conditions F(0) = 1 and Y'(Q) = 0. 
Equation l|ll[l assumes C > 0: a supercritical bifurcation. 
We have found that, at fixed 77, C > on an interval 
f-(rj) < f < f+(rf) that lies within the spinodal interval 
(hih)- On the intervals fi < f < f- and /+ < / < 
f2 the coefficient C becomes negative which indicates a 
subcritical bifurcation. The solid line in Fig. 6 shows the 
supercritical bifurcation curve l|ll|l for 77 = 11,050 and 
/ = 0.025. Here A c ~ 0.514 and T ~ 0.142. 
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When A is well above A c , the weakly nonlinear theory 
is invalid, and a numerical solution of the fully nonlinear 
hydrostatic problem (0}-© is needed for the determina- 
tion of \Y C \. An alternative approach is a hydrodynamic 
simulation, that is a numerical solution of the hydrody- 
namic equations. Numerical simulations of this type were 
done in Ref. ^(| for a different version of constitutive 
relations |10| and a different boundary condition at the 
driving wall. It was observed that the phase-separation 
instability produces multiple clusters whose further dy- 
namics proceed as gas-mediated competition and coars- 
ening. Direct merging of clusters can also occur. The 
final symmetry-broken state, as observed in the hydro- 
dynamic simulations, is always a single, almost densely 
packed stationary 2D cluster coexisting with gas (or di- 
lute bubble coexisting with denser fluid). The cluster 
is located in one of the system's corners (unless peri- 
odic boundary conditions are used). This scenario was 
confirmed in a hydrodynamic simulation of the present 
system (for r? = 11,050, / = 0.025 and A = 3) done by 
E. Livne j^. A density map of the hydrodynamic final 
state in this case is shown in Fig. 2D. The steady-state 
value \Y C \ ~ 0.265, obtained in this simulation, is shown 
by the circle in Fig. 6. 



IV. EMD SIMULATIONS 

A. Simulation method, parameters and diagnostics 

We put the predictions of the granular hydrostatics 
into test by doing extensive EMD simulations of this sys- 
tem. Most of the simulations were done with N = 2 - 10 4 
particles: hard disks of diameter d = 1 and mass m = 1. 
The thermal wall temperature is To = 1, so the scaled 
time unit is d (ro/To) 1 / 2 = 1. A standard event-driven 
algorithm |24| was used. Two of the hydrodynamic pa- 
rameters, n = 11, 050 and / = 0.025, were fixed in all sim- 
ulations, while A was varied in the range of 0.1 < A < 3. 
This was achieved by varying L x , L y and r. Indeed, for a 
fixed 77, /, A and N the coefficient of normal restitution, 



r = 1 - 



and the system's dimensions, 
1/2 



V3N 
~2fA 



1/2 



and L y = 



(12) 



(13) 



are uniquely determined. For the values of the param- 
eters that we used, r was always in the range of nearly 
elastic collisions: r > 0.977. The initial spatial distri- 
bution of the particles was (statistically) uniform, while 
the initial velocity distribution was Maxwell's with the 
wall temperature Tq = 1. The center-of-mass coordi- 
nate Y c (t) was used as a quantitative probe of the lateral 
asymmetry of the system. Before taking the steady-state 



measurements we waited until transients died out. This 
was monitored by the time-dependence of the average 
kinetic energy of the particles (that first decayed and 
then approached an almost constant value) and by the 
time-dependence of the center of mass itself, see below. 
Selected movies of these simulations can be downloaded 



from http : / /bioinf . charite . de/kies/ giantfluctuations / 



B. Final states at different A 

The EMD simulations showed that, at aspect ratios 
well below the threshold value of A = A c ~ 0.512, the 
final state is a (weakly fluctuating) stripe state. The 
number density profile versus x, found in the simula- 
tions, compares very well with the hydrostatic solution 
(see Fig. 1), while Y c {t) stays close to zero. Notice that 
the Jenkins-Richman constitutive relations [2(j, that we 
used in this comparison, do not include any fitting pa- 
rameters. Therefore, well below the instability threshold 
in A, the hydrostatic solution yields a quantitatively ac- 
curate leading-order description of the system. 

At aspect ratios well above the instability threshold we 
always observed several clusters nucleating at the wall op- 
posite to the driving wall. The cluster dynamics (Fig. 2 
A to C) proceeds as gas-mediated competition and coars- 
ening (sometimes as direct mergers) of clusters, in accord 
with hydrodynamic simulations [la ]. As time increases, 
the number of clusters goes down, and only one dense 
cluster, fluctuating around its average position in one of 
the two corners, opposite to the thermal wall, finally sur- 
vives. Fig. 2C shows a snapshot of the final state for 
A = 3. For comparison, Fig. 2D shows a density map 
of the final steady state obtained by E. Livne in a hy- 
drodynamic simulation for the same hydrodynamic pa- 
rameters. The center-of-mass position Y c of the steady 
state agrees well with the average-in-time center-of-mass 
position, measured in the EMD simulations, as shown 
by the circle in Fig. 6. This indicates that, well above 
the instability threshold, the hydrostatic theory describes 
the steady states of the system well. We can also refer 
the reader to the recent EMD simulation results for very 
large aspect ratios 17] . As no appreciable fluctuations 
around a broken-symmetry steady state were reported, 
one can safely assume that the broken-symmetry steady 
states observed in Ref. 01 should be also describable by 
the hydrostatic theory. 

The system behavior changes dramatically, however, as 
the aspect ratio A approaches A c . We found that, in a 
wide region of A around A c , the final state of the system 
exhibits large-amplitude irregular oscillations, as dense 
clusters at the wall opposite to the driving wall nucle- 
ate, move in the lateral direction, dissolve and reappear. 
Figure 3 shows a typical sequence of snapshots from an 
EMD simulation for A = 1. 

Figure 4 shows the time history of the center-of-mass 
coordinate Y c for six different values of A. One can see 
that, in a wide region of intermediate A, the center- 
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of- mass coordinate Y c (t) shows large- amplitude irregu- 
lar oscillations. Noticeable are multiple zero crossings of 
Y c {t) at aspect ratios above the hydrodynamic bifurcation 
point A c ~ 0.512 (Fig. 4 c-e). Smaller but still signifi- 
cant irregular oscillations are also observed below A c , as 
if the system persistently tends to break the lateral sym- 
metry there. The hydrostatic picture is recovered when 
one moves farther away, in any direction, from the region 
of A ~ A c . Indeed, Fig. 4e shows that zero crossings 
of Y c (t) occur less often for A = 1.3, than for A = 0.7 

A B C D 




or 1. At still larger A (Fig. 4f) no zero crossings are 
observed for any reasonable simulation time, and Y c fluc- 
tuates around a constant value that is very close to that 
predicted by the hydrostatic theory (and shown by the 
circle in Fig. 6). 

To better characterize the fluctuation-dominated re- 
gion, we computed the probability distribution function 
P(|Y C |) of different values of \Y C \ in a statistical steady 
state, that is, after transients die out. The stationarity 
of the remaining data was tested by dividing the respec- 
tive time interval into three sub-intervals and checking 
that the differences in P(|Y C |) for the sub-intervals are 
small and not systematic. The probability distribution 
P(|F C |) is shown, at different A, in Fig. 5. At A < A c 
the maximum of P(|Y C |) is at \Y C \ — 0, and it is rela- 
tively narrow. Correspondingly, there is no symmetry 
breaking there, the fluctuations are relatively small, and 
the hydrostatic theory yields an accurate leading-order 
description. At A 3> A c , the maximum of P(|Y C |) is at a 
non-zero \Y C \. This is a clear manifestation of symmetry- 
breaking: a dense cluster develops in one of the corners 
away from the driving wall. The probability distribution 
P(|Y C |) is also quite narrow here, the fluctuations are rel- 
atively small, and there is a good agreement between the 
hydrostatic theory and EMD-simulations. On the con- 
trary, in a wide region of A around A c , the probability 
distribution P(|Y C |) is very broad, and the hydrostatic 
theory breaks down. By following the position of the 



FIG. 2: Nucleation and coarsening of clusters as observed in 
an EMD simulation with N = 2 • 10 4 particles for Tj = 11, 050, 
/ = 0.025 and A = 3. The hot wall is on the right. The 
scaled times are 14,425 (A), 26,218 (B) and 191,616 (C). 
Figure D is a density map of the steady state obtained by E. 
Livne in a simplified hydrodynamic simulation for the same 
hydrodynamic parameters 
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FIG. 3: Irregular lateral cluster dynamics for A = 1 as ob- 
served in an EMD simulation with N = 2 • 10 4 particles for 
r\ — 11,050 and / = 0.025. The time progresses from left to 
right, starting from the upper row. The hot wall is on the 
right. 
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FIG. 4: Y c versus time for 77 = 11,050 and / = 0.025 and 
different values of the aspect ratio A, as observed in EMD 
simulations with N = 2 • 10 4 particles. Time here is propor- 
tional to the number of particle collisions; t = 500 corresponds 
to 505, 036 scaled time units. 



maximum of P(|Y C |) at different A (see Fig. 6), one can 
see that the symmetry-breaking transition occurs some- 
where in the region of 0.3 < A < f.0. Because of the 
extreme flatness and broadness of the probability distri- 
bution P(|Y^|) in this region, a more accurate estimate of 
the position of the maximum of P(|Y C |) requires a much 
better statistics (that is, a much longer simulation time) 
than we could afford in this series of simulations [26( . 
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FIG. 6: Fhe effective bifurcation diagram of the system for 
77 = 11,050 and / = 0.025, observed in EMD simulations 
with with N — 2 ■ 10 4 particles. Diamonds show, for each A, 
the positions of the maxima of the probability distribution 
function P(|Y C |). Above the transition, the error bars show 
the errors in the estimation of the position of the maximum of 
P(|y c |). Below the transition the error bars show the errors in 
the estimation of (Y c ) : the time average of Y c . The solid line is 
the bifurcation curve (II II) close to threshold. The empty circle 
at A = 3 shows the result of the hydrodynamic simulation by 
E. Livne. The dashed line is an interpolation between the 
solid line and the empty circle. 



FIG. 5: The probability distribution function P(|Y C |) of the 
fluctuating final state of the system for rj = 11, 050 and / = 
0.025 and different values of the aspect ratio A, as observed in 
EMD simulations with N = 2 • 10 4 particles. In order to show 
all the graphs on the same scale, the probabilities (rather than 
the probability densities) for each bin are shown. 

Noticeable in Fig. 6 is a systematic discrepancy, within 
the wide fluctuation-dominated region, between the po- 
sitions of the maxima of P(|Y C |) and the hydrostatic bi- 
furcation curve computed in Sec. III. We even cannot ex- 
clude a change in the character of bifurcation caused by 
the fluctuations (apparently without shifting the bifurca- 
tion point). Indeed, the maxima of P(|Y C |) at A = 1.0, 
1.3 and 2.0 appear to lie on a straight line passing through 
the theoretical transition point A c ~ 0.5. As A increases 
further, the discrepancy between the positions of the 
maxima of P(| y r |) and the theoretical bifurcation curve 
goes down [25J . Importantly, the fluctuation-dominated 
region 0.3 < A < 1.0 does include the hydrostatic tran- 
sition point A c ~ 0.5. 

We should stress that the failure of hydrostatics is ob- 
served at intermediate values of the aspect ratio A, when 
the hydrodynamic parameters 77 and /, and the number 
of particles N, are fixed. In view of Eq. (|12|) . while in- 
creasing A, one increases the inelasticity of particle col- 
lisions 1 — r. That the hydrostatic theory fails at inter- 
mediate values of the inelasticity, and improves at small 
enough, or large enough inelasticities, excludes the in- 
elasticity itself as the reason for the failure. 



C. Simulations with different N 

We did a series of simulations with different number 
of particles N in order to verify the hydrodynamic scal- 
ing and investigate the A"-dependence of the (relatively 
weak) fluctuations well below and well above A c . These 
additional simulations were done for A = 0.1 and three 
values of N: 5 • 10 3 , 10 4 and 1.5 • 10 4 , and for A = 3.0 
and AT = 4 • 10 4 . 

When varying A" at fixed A, we kept the hydrody- 
namic parameters r\ = 11,050 and / = 0.025 constant. 
Therefore, if the hydrostatic equations provide a correct 
leading-order theory of the steady states far below and 
far above A c , the time-averaged steady state values of 
Y c should become A"-independent for large enough N. 
Figure 7 shows Y c versus time for A = 0.1 at the four 
different values of N. One can see that, in all these cases, 
the average value of Y c is close to zero as expected, while 
fluctuations are relatively small. Figure 8 shows the dy- 
namics of Y c (t) for A = 3 and two different values of N: 
2- 10 4 and 4- 10 4 . Here the symmetry-breaking is evident, 
as a dense cluster develops in a corner. With a moderate 
accuracy determined by the relatively high level of fluc- 
tuations of Y c , the average values of Y c at late times are 
close to each other. Therefore, well below and well above 
A c the hydrodynamic scaling is obeyed. 

Simulations with fixed scaled parameters 77, / and A, 
but different N can also help in identifying the mecha- 
nism of breakdown of the hydrostatic theory at aspect 
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FIG. 7: Y c versus time for 77 = 11, 050, / = 0.025 and A = 0.1, 
for N = 5000 (a), 10 000 (b), 15 000 (c) and 20 000 (d), as 
observed in EMD simulations. Time units are the same as in 
Fig. 4. 



ratios around A c . Indeed, it is natural to interpret the 
giant oscillations, shown in Fig. 4c-e, in terms of a strong 
coupling between the two bifurcated states predicted by 
the hydrostatic theory. One possible scenario of this cou- 
pling (which we call Scenario I) relies on the discrete- 
particle noise, unaccounted for by granular hydrodynam- 
ics. Below A c , the discrete-particle noise is expected to 
cause fluctuations, that is to broaden the distribution of 
Y c as indeed observed in Fig. 5. If Scenario I is correct, 
the standard deviation a of Y c (t) from its average value 
should vanish as N goes to infinity, at fixed hydrody- 
namic parameters 77, / and A. 

Another possibility (Scenario II) is that the fluctua- 
tions persist in the limit of N — > 00. If this is the 
case, the dominating mechanism of fluctuations has a 
purely hydrodynamic nature and should be explainable 
by a full hydrodynamic analysis (as opposed to our hy- 
drostatic analysis, and to the simplified hydrodynamic 
simulations that used a model Stokes friction instead of 
the full viscosity). Here the coupling between the two 
symmetry-broken states may be due to either an unstable 
hydrodynamic mode (Scenario Ha), or a weakly damped 
mode (Scenario lib). In Scenario lib, a should vanish, as 
N — > 00, if one waits for a sufficiently long time. There- 
fore, to distinguish between the two sub-scenarios, one 
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FIG. 8: Y c versus time for 77 = 11, 050, / = 0.025 and A = 3, 
for two different values of N, as observed in EMD simula- 
tions. The thick line corresponds to N = 2 ■ 10 4 , the thin line 



corresponds to N 
Fig. 4. 



4 • 10*. Time units are the same as in 



should, in addition to the limit of N — > 00, take the limit 
of t —>■ 00. 

Obviously, one is unable to take any of these two limits 
in actual EMD simulations, where the maximum achiev- 
able values of N and t are limited by the available com- 
puter resources. So what was observed in our EMD sim- 
ulations with different TV? Figures 7 and 9 show what 
happens well below A c , when N increases from 5 000 
to 20 000. One can see from Fig. 7 that, as N grows, 
the high-frequency components of the fluctuations do de- 
crease, but the low frequency component does not show 
any pronounced decrease. Overall, the fluctuation spec- 
trum moves towards the lower frequencies. As the result, 
a good resolution of the low-frequency part of the power 
spectrum requires longer and longer simulations (which 
rapidly become prohibitively long). This introduces an 
additional, non-trivial constraint on simulations with a 
large number of particles. A similar situation occurs well 
above A c . Figure 8 does indicate that a goes down as 
N goes up from 20 000 to 40 000. However, one also ob- 
serves that, as N grows, the role of the low- frequency 
components of the fluctuations increases. 

Hydrodynamics provides a hint for the mechanism of 
the "red shift" of the power spectrum with an increase 
of N. There are four hydrodynamic modes in the sys- 
tem: two acoustic modes, the entropy mode and the 
shear mode. The frequencies of the acoustic modes are 
the highest, as they are determined by the "ideal" (non- 
dissipative) terms in the hydrodynamic equations, and 
they scale like the inverse system size. The frequencies of 
the entropy and shear modes are much lower, as they are 
determined by the transport coefficients: the heat con- 
duction, viscosity and inelastic loss rate, and they scale 
like the inverse square of the system size. In the units of 
d = m = Tq = 1, and at fixed hydrodynamic parameters 
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77, / and A, a larger A implies a larger system, see Eqs. 
(|13fl . Correspondingly, as A increases, the characteristic 
frequencies of the entropy /shear modes go down much 
faster than those of the acoustic modes. Therefore, it 
seems likely that one of these modes is responsible for 
the low- frequency components of the fluctuations. A re- 
lated issue is that, in contrast to the hydrostatic problem 
|JT)|. the full time-dependent hydrodynamic problem has 
an additional scaled parameter: d/L x . This parameter 
describes the role of the dissipative terms compared to 
the "ideal" terms in the hydrodynamic equations. As it 
is clear from Eq. H13|) . when increasing A at constant 77 
and /, one reduces this additional parameter. There- 
fore, as A increases, the low- frequency shear/entropy 
modes should become more and more persistent. As 
these modes are not necessarily broad-band, a might 
cease to provide a good characterization of the system 
at large A. 

Still, if one continues following a as A increases, one 
observes (see Fig. 9) that a decreases much slower than 
the classic dependence N~ x l 2 characteristic of equilib- 
rium systems. If one attempts to interpret the decrease 
of a with an increase of A in terms of an empiric power 
law, one obtains an exponent —0.23, instead of the classi- 
cal value of —1/2 for equilibrium systems. Importantly, 
we did reproduce the classical A -1 / 2 scaling of a in a 
control series of simulations with the same / and A, but 
with 77 = (elastic collisions). Moreover, a good quanti- 
tative agreement was obtained with a theoretical result 
for er that directly follows from the classic expression for 
the density correlation function in equilibrium [2^. We 
also found that, for the same total number of particles 
A, the fluctuation levels in the elastic case are signif- 
icantly lower than in the inelastic case. That is, well 
below A c , the fluctuations, though much smaller than 
those observed for A ~ A c , are still large compared to 
the elastic case. 

To summarize this subsection, our simulations with dif- 
ferent A strongly indicate that the hydrostatic equations 
provide a correct leading order theory of this system well 
below and well above A c . On the other hand, the simula- 
tions proved to be insufficient for determining the mecha- 
nism of giant fluctuations that we observed in this system 
at A ~ A c . We cannot even be sure at this point whether 
the fluctuations (or, more precisely, their low-frequency 
components) persist or not as A — ► 00. 



V. SUMMARY AND DISCUSSION 

The main results of this work can be summarized in 
the following way. Granular hydrostatics, in combination 
with simplified hydrodynamic simulations, correctly pre- 
dict the phase separation instability in this prototypical 
driven granular system. Well above and well below the 
critical value of the aspect ratio A c , the hydrostatic the- 
ory describes the steady state of the system well. How- 
ever, in a wide region of aspect ratios around A c the 
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FIG. 9: Shown is a, the standard deviation of Y c from its 
(almost zero) average value, versus the number of particles N 
for 77 = 11,050, / = 0.025 and A = 0.1. The symbols show 
the simulation results. The curve shows, as a reference, the 
power-law dependence a = BN -13 with exponent — —0.23, 
see the text. 



system is dominated by fluctuations, and the hydrostatic 
theory fails. The fluctuation levels are anomalously high 
even relatively far from the hydrostatic bifurcation point, 
and they certainly do not exhibit the classic A -1 / 2 scal- 
ing with the number of particles A. 

Though we are unable to pinpoint the mechanism of 
excitation of the giant fluctuations, we can suggest two 
different scenarios for their origin. In Scenario I the fluc- 
tuations are driven by discrete particle noise. Indeed, 
it is well known that discrete particle noise can drive 
relatively large fluctuations in the vicinity of thresholds 
of hydrodynamic instabilities |T^ and non-equilibrium 
phase transitions [2^. Fluctuations of this type should 
vanish as one increases indefinitely the number of parti- 
cles in the system, keeping the hydrodynamic parameters 
constant. Unfortunately, our simulations with different 
A, but fixed 77, / and A, have been insufficient to prove 
or disprove this scenario. 

A difficulty with Scenario I is that the fluctuations are 
so big in so wide a region of aspect ratios. No anomaly 
of this type has been observed in any other symmetry- 
breaking instability of granular flow, even with much 
smaller numbers of particles. As an example, let us con- 
sider for a moment the same system, but introduce grav- 
ity in the x direction. Now the granular gas is heated 
from below, and the system exhibits another symmetry- 
breaking instability: thermal convection, similar to the 
Rayleigh-Benard convection of classical fluids. The tran- 
sition to convection occurs via a supercritical bifurcation 
Though EMD simulations of thermal gran- 
ular convection [23 involved only A = 2, 300 particles 
(which is much less than A = 2 • 10 4 used in the present 
work) , a sharp supercritical bifurcation was observed, in 
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agreement with a hydrodynamic analysis [30l l3l| . By 
comparison, the giant fluctuations, observed in a wide 
region of A in the present work, are an anomaly, as one 
needs some (hydrodynamic?) mechanism of strong am- 
plification of the discrete-particle noise. 

If Scenario I proves to be correct, the corresponding 
theory can be developed in the framework of Fluctuat- 
ing Hydrodynamics |27| . generalized to granular gases 
in the limit of nearly elastic collsions. Fluctuating Hy- 
drodynamics is a Langevin-type theory that takes into 
account the discrete character of particles by adding 
delta-correlated noise terms in the momentum and en- 
ergy equations j2]j. Fluctuating Hydrodynamics is by 
now well established for classical fluids in 3D, including 
non-equilibrium states 0, • We should mention here 
that the 2D case has an additional difficulty The cou- 
pling of fluctuations here is anomalously strong, even in 
the elastic case: the transport coefficients diverge in the 
thermodynamic limit, except for a sufficiently dilute gas 
[33j. Therefore, one can hope to generalize the Fluc- 
tuating Hydrodynamics to the 2D gas of inelastic hard 
spheres in the dilute limit [34|]. Close to the phase sepa- 
ration threshold, the dilute limit holds with a reasonable 
accuracy. It would be interesting to investigate the phase 
separation problem in 3D, where important differences in 



the fluctuation behavior may occur. 

Alternatively, in Scenario II the low-frequency compo- 
nent of the giant fluctuations has a purely hydrodynamic 
origin and is driven either by a presently unknown hy- 
drodynamic instability (Scenario Ha), or by a long-lived 
transient mode (Scenario lib). Effects of these type are 
obviously missed by a hydrostatic analysis. They may 
have also been missed by the time-dependent hydrody- 
namic simulation p3| that employed a model Stokes fric- 
tion, rather than the hard-sphere viscosity, to acceler- 
ate the convergence to a steady state. If Scenario II is 
correct, the low- frequency component of the fluctuations 
should be observable in hydrodynamic simulations with 
the true hard-sphere viscosity. These simulations, there- 
fore, should be an important next step in the analysis of 
this fascinating problem. 
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